Stimulation experiments

Analysis of stimulation experiments.

library(Seurat)
library(dplyr)
library(ggplot2)
library(scClustViz)

Loading required package: shiny

Load integrated object. This object has been integrated with CCA, because there was lots of batch effect between the 0hr tissue and the other two experiments (which merged quite well).

load("~/Dropbox/Zoe/scf_version/analysis/healthy_sc/seurat_objects/no_dropletQC/3391/integrated_3391_stimulation_pmaio_unstim_downsampled_cca_kanchor5_scClustViz.RData")
# Set identity
Idents(scSeurat) <- "integrated_snn_res.0.4"

Check out the UMAP of clusters, original identities, and refined SCINA labels.

DimPlot(scSeurat, label = TRUE)

DimPlot(scSeurat, group.by = "tissue_type", label = FALSE)

DimPlot(scSeurat, group.by = "scina_labels", label = TRUE) & NoLegend()

pdf(file = './figures/stimulation/downsampled_UMAP.pdf')
DimPlot(scSeurat, label = TRUE)
dev.off()

null device 
          1 

pdf(file = './figures/stimulation/downsampled_tissue_type_UMAP.pdf')

DimPlot(scSeurat, group.by = "tissue_type") +
  ggtitle(NULL)
dev.off()

null device 
          1 

Look at tissue-type per cluster on a proportional level:

pt <- table(Idents(scSeurat), scSeurat$tissue_type)
pt <- as.data.frame(pt)
plot <- ggplot(pt, aes(x = Var1, y = Freq, fill = Var2)) +
  theme_bw(base_size = 15) +
  geom_col(position = "fill", width = 0.5) +
  xlab("Sample") +
  ylab("Proportion") +
  theme(legend.title = element_blank())
plot
pdf(file = './figures/stimulation/downsampled_tissue_type_barplot.pdf')

plot
dev.off()

png 
  2 

Look for T-cell signatures.

FeaturePlot(scSeurat, features = c("CD3E", "NKG7"))

FeaturePlot(scSeurat, features = c("LEF1", "XCL1;XCL2"))

FeaturePlot(scSeurat, features = c("TOX", "CD8A"))

FeaturePlot(scSeurat, features = c("CCR7"))

Look for macrophage signatures:

FeaturePlot(scSeurat, features = c("LGALS1", "TYROBP"))

FeaturePlot(scSeurat, features = c("MARCO", "C1QC"))

FeaturePlot(scSeurat, features = c("LYZ-1", "S100A8"))

FeaturePlot(scSeurat, features = c("DEFA6;DEFA4;DEFA5", "CTSS"))

FeaturePlot(scSeurat, features = c("FLT3", "S100A12"))

Look for B and plasma cell signatures

FeaturePlot(scSeurat, features = c("CD79B", "JCHAIN"))

Cluster 12 is the only cluster that is obviously dominated by the unstimulated experiments. What cells are these? Looks like they are hepatocytes. A lot of the markers are periportal, but all hepatocytes are in this cluster.

# Isolate DE genes
clust11 <- sCVdata_list$res.0.4@DEvsRest$`12`[order(sCVdata_list$res.0.4@DEvsRest$`12`$Wstat,
                                                    decreasing = TRUE), ]
rownames(clust11)[1:50]

 [1] "PAH"                                      "LOC114082496"                             "TDO2"                                     "MAT1A"                                    "GYS2"                                    
 [6] "SLCO1B3;SLCO1B3-SLCO1B7;SLCO1B7;SLCO1B1"  "FBP1"                                     "CPS1"                                     "SLC2A2"                                   "ALDOB"                                   
[11] "MTSS1"                                    "IGFBP1"                                   "CYP3A4;CYP3A7-CYP3A51P;CYP3A7-2"          "mikado.WCK01-AAH20201022-F8-rsc00143G668" "PCK1"                                    
[16] "SHMT1"                                    "CYP3A4;CYP3A7-CYP3A51P;CYP3A7"            "ARG1"                                     "CYP3A4;CYP3A7-CYP3A51P;CYP3A7-3"          "LOC107139912"                            
[21] "ALB-1"                                    "SLC38A4"                                  "LOC107139914"                             "SORBS2"                                   "TAT"                                     
[26] "C5-1"                                     "KLKB1"                                    "Slc22a19-1"                               "G6PC1"                                    "RBP4"                                    
[31] "AGMO"                                     "HSD17B11"                                 "CMBL"                                     "SLC25A13"                                 "ZBTB16"                                  
[36] "C8G"                                      "PYGL"                                     "HSD11B1"                                  "mikado.WCK01-AAH20201022-F8-rsc00060G101" "CES3"                                    
[41] "MTARC2"                                   "IDO2"                                     "AQP9"                                     "ACACB"                                    "GOT1"                                    
[46] "CFI"                                      "NR1I2"                                    "ABCB11"                                   "SOX5"                                     "BHMT2"                                   

FeaturePlot(scSeurat, features = c("PAH","PCK1"))

FeaturePlot(scSeurat, features = c("CYP2E1","FETUB"))

FeaturePlot(scSeurat, features = c("ACACB","ELOVL6"))

FeaturePlot(scSeurat, features = c("POLR2D", "ALB-1"))

Let’s look at some of the key markers that Sonya sent: IL2, IL4, IL6, IL7, IL10, IL12A, TGFB

FeaturePlot(scSeurat, features = c("IL2", "IL4"))

FeaturePlot(scSeurat, features = c("IL6", "IL6-1"))

FeaturePlot(scSeurat, features = c("IL7", "IL10"))

FeaturePlot(scSeurat, features = c("IL12A", "TGFB1"))

It’s pretty hard to see which cytokines are expressed in which samples, so I’m going to try splitting up the samples.

DimPlot(scSeurat, split.by = "tissue_type", label = TRUE) & NoLegend()
pdf(file = './figures/stimulation/downsampled_UMAP_split.pdf')

DimPlot(scSeurat, split.by = "tissue_type", label = TRUE) & NoLegend()
dev.off()

png 
  2 

Now let’s take a look at the features on this split plot.

genes <- c("IL2","IL4","IL6","IL7","IL10","IL12A","TGFB1","TGFB2","TGFB3","TNF","IL18","IL2RA","CCL5","IFNG","LIF")
for (gene in genes) {
  print(FeaturePlot(scSeurat, feature = gene, split.by = "tissue_type"))
  pdf(paste("figures/stimulation/downsampled_UMAP_split_", gene, ".pdf", sep = ""))
  print(FeaturePlot(scSeurat, feature = gene, split.by = "tissue_type"))
  dev.off()
}

Warning: All cells have the same value (0) of “IL2”

More plots

Cluster 2 looks like it generally has more IL genes expressed. Isolate DE genes for this cluster.

clust1 <- sCVdata_list$res.0.4@DEvsRest$`2`[order(sCVdata_list$res.0.4@DEvsRest$`2`$Wstat,
                                                    decreasing = TRUE), ]
cat(rownames(clust1)[1:50])

PTPRC TOX PRKCH INPP4B mikado.WCK01-AAH20201022-F8-rsc00190G443 ARID5B IKZF3 ARHGAP15 SATB1 DOCK10 TRPS1 RASA2 SAMSN1 NSD3 MBNL1 SH3KBP1 EMB SYTL3 ICOS mikado.WCK01-AAH20201022-F8-rsc00083G981 mikado.WCK01-AAH20201022-F8-rsc00275G164 THEMIS MSI2 SLAMF1 PITPNC1 NXF2;NXF2B;NXF5 SRSF7 IKZF1 RAPGEF6 STAT4 SLC38A1 RAP1A RABGAP1L ZC3HAV1 BCL2 TESPA1 APBB1IP GIMAP7 DHX8 DDX5 UBAC2 HUWE1 CD3E CTLA4 CD44 GIMAP1 RGS1 MTHFD1L RIPOR2 TANK

Visualise these genes in violin plots.

Interpretation: Now I see that IL2 and IL4 are expressed in more than just cluster 2, but it is nearly always in PMAIO. Cluster 8 is definitely where IL6 is expressed. IL7 is a bit more spread out across treatment types, but it’s definitely present in PMAIO in clusters 4, 5, and 11. Cluster 11 is plasma cells, but I’m honestly not sure of 4 and 5 - low signal hepatocytes? IL10 is in clusters 2 and 11 PMAIO. IL12A is a bit more in 4 and 5 PMAIO.

VlnPlot(scSeurat, feature = "IL2", split.by = "tissue_type") # Only in PMAIO cluster 2

VlnPlot(scSeurat, feature = "IL4", split.by = "tissue_type") # Only in PMAIO cluster 2

VlnPlot(scSeurat, feature = "IL6", split.by = "tissue_type") # Stronger in endo and mesem PMAIO

VlnPlot(scSeurat, feature = "IL6-1", split.by = "tissue_type") # No signal

VlnPlot(scSeurat, feature = "IL7", split.by = "tissue_type") # More in PMAIO cluster 4/5

VlnPlot(scSeurat, feature = "IL10", split.by = "tissue_type") # More in PMAIO cluster 2

VlnPlot(scSeurat, feature = "IL12A", split.by = "tissue_type") # A bit more in 1,4,5

VlnPlot(scSeurat, feature = "TGFB1", split.by = "tissue_type") # Pretty even across all

VlnPlot(scSeurat, feature = "TNF", split.by = "tissue_type")

Make Violin plots with ggplot

genes <- c("IL2","IL4","IL6","IL7","IL10","IL12A","TGFB1","TGFB2","TGFB3","TNF","IL18","IL2RA","CCL5","IFNG","LIF")
for(gene in genes) {
  vln_df = data.frame(gene = scSeurat[["RNA"]]@data[gene,], 
                      Cluster = Idents(scSeurat), 
                      tissue_type = scSeurat$tissue_type)
  plot <- ggplot(vln_df, aes(x = Cluster, y = gene)) + 
    geom_violin(aes(fill = tissue_type), trim=TRUE, scale = "width") +
    ylab(gene) +
    theme_bw()
  print(plot)
  pdf(file = paste("figures/stimulation/downsampled_vlnPlot_split_", gene, ".pdf", sep = ""))
  print(plot) # Need print function because R doesn't autoprint in loops
  dev.off()
}

Looking for more genes

gene <- "IL18"
vln_df = data.frame(gene = scSeurat[["RNA"]]@data[gene,], 
                      Cluster = Idents(scSeurat), 
                      tissue_type = scSeurat$tissue_type)
  plot <- ggplot(vln_df, aes(x = Cluster, y = gene)) + 
    geom_violin(aes(fill = tissue_type), trim=TRUE, scale = "width") +
    ylab(gene) +
    theme_bw()
  print(plot)

Perform DE test on T-cells from PMAIO vs unstim T-cells. I see TNF, IL2RA, but I also see hepatocyte genes like ALB-1, STAT3

cat(row.names(de_results)[1:100], sep = " ")

LOC107134468 CCL4L2;CCL4L1;CCL4 mikado.WCK01-AAH20201022-F8-rsc00073G191 MTHFD1L NFKB1 NXF2;NXF2B;NXF5 CCND3 IFNG SATB1 SKAP1 ARID5B mikado.WCK01-AAH20201022-F8-rsc00334G1 CCL3L1;CCL3L3;CCL3;CCL18 SLAMF1 CD200 GIMAP7 LOC114086011-6 GBP7;GBP4 CBLB TRAF3IP3 NR4A3 ALB-1 TNFRSF9 SLC25A26 SYNE1 STAT3 UTRN HUWE1 PDE7A IL2RA PTPRCAP APOA1 TNRC6B Ifi214;Ifi209 ETS1 ANKRD44 MAPK6 PDE3B RESF1 CD27 mikado.WCK01-AAH20201022-F8-rsc00138G153 PICALM JARID2 MDN1 CHIC1 CCL3L1;CCL3L3;CCL3;CCL18-2 CKLF;CKLF-CMTM1 HECA FOXN3 MED13L SEC23B ZEB2 AMN1 CD3G SEC14L1 BATF LNPEP IRF4 ALB LOC107134769-1 IQGAP2 LOC114086011-12 SERPINA1 mikado.WCK01-AAH20201022-F8-rsc00234G885 TMEM156 STK4 ATP8A1 PECAM1 FAM13A MTHFD2 CELF2 SAMSN1 NUCB2 GIMAP1 SND1 VMP1 MLLT3 ASAP1 RBP4 TNF CRTAM FYN ACTB mikado.WCK01-AAH20201022-F8-rsc00190G443 PPP3CC HBS1L TTC27 B4GALT5 DDX5 ARHGDIB LAT CXCL2;CXCL3;CXCL1-1 ITM2B SENP7 DOCK8 RIPOR2 APBB1IP TRAF3 SAMHD1 CD99